#####
#####
###    Create Figure 1 in supplemental report, 
###    "Assessing Multivariate Normality in
###    Bootstrap Replications of the MDS
###    Point Coordinates." 
#####
#####

library(lattice)


###
###   Read bootstrap replication data
###   from file, "Bootstrap estims, 
###   with names.txt"
###

bsreps <- read.table(file.choose(), header = T)


#####
#####
###     Create trellis display containing a scatterplot
###     of the sorted D-squared values for each candidate
###     versus the quantiles of the Chi-squared(2)
###     distribution.
#####
#####


###
###   Create data frame with name and dimensions only
###

cands <- bsreps[, c("name", "dim1", "dim2")]

###
###   Create object with sorted quantiles of 
###   chi-squared(2) distribution
###

quant.chisq2 <- qchisq((1:50 - .5)/50, df = 2)

###
###   Create vector of names
###

cnames <- levels(cands$name)

full.cnames <- c("Ashcroft", "Bill Clinton", "Cheney", 
   "Democratic Pty", "Edwards", "George Bush", 
   "Hillary Clinton", "Kerry", "Laura Bush", "McCain",
   "Nader", "Powell", "Republican Pty")
   

   for (i in 1:length(cnames)) {
      curr.cand <- cands[cands$name == cnames[i], 2:3]   
      curr.means <- mean(curr.cand)
      curr.covar <- cov(curr.cand)
      curr.centered <- as.matrix(t(apply(curr.cand, MARGIN = 1, function (x) {x - curr.means})))
      curr.quants <- 
         apply(curr.centered, MARGIN = 1, 
         function (x) {t(as.matrix(x)) %*% solve(curr.covar) %*% as.matrix(x)})
      if (i == 1) {cand.mnorm <- data.frame(name = full.cnames[i], quant.cand = sort(curr.quants), quant.chisq2)}
      if (i > 1) {cand.mnorm <- rbind(cand.mnorm, data.frame(name = full.cnames[i], quant.cand = sort(curr.quants), quant.chisq2))}
      }

###
###   The following function creates Figure 1
###   in the supplemental report. Note that it
###   removes a single outlier with the "subset=" argument
###

      
xyplot(quant.cand ~ quant.chisq2 | name, data = cand.mnorm,
   aspect = 1,
   panel = function (x, y) {
   panel.abline(a = 0, b = 1, col = "black")
   panel.xyplot(x, y, cex = .75)
   },
   subset = quant.cand < 18,
   par.strip.text = list(cex = .75),
   scales = list(at = c(1, 6, 11), tck = .6, cex = .75),
   xlim = c(-.75, 12.75),
   ylim = c(-.75, 12.75),
   xlab = expression(chi[(2)]^2~~"quantiles"),
   ylab = expression("Mahalanobis"~~italic(D)^2~~"for bootstrap replications")
)

